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Abstract. We consider several solitons moving in a slowly varying external field. We 
show that the effective dynamics obtained by restricting the full Hamiltonian to the finite 
dimensional manifold of Af-solitons (constructed when no external field is present) provides 
a remarkably good approximation to the actual soliton dynamics. That is quantified as 
an error of size h 2 where h is the parameter describing the slowly varying nature of the 
potential. This also indicates that previous mathematical results of Holmer-Zworski [8] 
for one soliton are optimal. For potentials with unstable equilibria the Ehrenrest time, 
log(l//i)//i, appears to be the natural limiting time for these effective dynamics. We 
also show that the results of Holmer-Perelman-Zworski [7 for two mKdV solitons apply 
numerically to a larger number of interacting solitons. We illustrate the results by applying 
the method with the external potentials used in Bose-Einstein soliton train experiments 
of Strecker et al [14]. 



In many situations a wave moving in a slowly varying field, that is, a field described 
by a potential whose derivatives are much smaller than the oscillations/ width of the wave, 
can be described using classical dynamics. This is the basis of the semiclassical/short wave 
approximation, perhaps best known in the case of the linear Schrodinger equation, 



where V is an infinitely differentiate potential. A typical result concerns a propagation of 
a coherent state 



maximally concentrated near the point (x ,Co) i n the position-momentum space. In that 

case, 



1. Introduction 



(i.i) 



ihd t u 





(1.2) 




< t < T(h) , 



where lmdl 



(p > 0, Im ip > 0, and 
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Exact solution; N = 4 



Effective Dynamics; N = 4 




- Exact Solution 

- Effective Dynamics 





Figure 1. A side- by-side comparison of the effective dynamics versus the 
exact solution of (1.6) for 4 solitons with the potential W(x) = — 100e cosx . 
The plot on the left shows the absolute value of the solutions up to time 
t — 1. The plot on the right shows the real part of the solutions at times 
t — 0.5 and t — 1. Compared to the solutions in Figures [4] and [5| much less 
discrepancy between the two solutions is visible. 



The time of the validity of (1.2), T(/i), depends on the properties of the flow (1.3), and in 



general it is limited by the Ehrenfest time. 

1 



(1.4) T{h) ~ log 



h 



see pQ for a recent discussion on the case of one dimension. 



The approximation (1.2) means that the solution is concentrated for logarithmically long 
times on classical trajectories. The phase tp and the amplitude a can be described very 
precisely and ao can be refined to give an asymptotic expansion - see [6] for an early 
mathematical treatment and [13] for more recent developments and references. 

In this paper we consider the Gross-Pitaevski equation, which is the cubic non-linear 
Schrodinger equation with a potential: 

(1.5) ihd t u = -^h 2 d 2 x u - u\u\ 2 + V(x)u . 

It provides a mean field approximation for the evolution of Bose-Einstein condensate in an 
external field given by the potential V(x) - see the monograph [12] and references given 
there. Questions about propagation of localized states are also natural in the setting of 



(1.5) and have been much studied. One direction is described in a recent monograph [2]. 



In this note we present a numerical study of multiple soliton propagation for (1.5) and 
show that it can be described very accurately using a natural effective dynamics - see 
Figure [T] That effective dynamics is based on mathematical results of Holmer-Zworski 
[8] and Holmer-Perelman-Zworski [7] and we refer to those papers for pointers to earlier 
mathematical works on that subject. 



EFFECTIVE DYNAMICS FOR iV-SOLITONS 



3 



Following the convention of earlier papers - see Frohlich et al [4J - we rescale equation 
( |1.5D so that the parameter h is in the potential which is now slowly varying: 

(1.6) id t u = — ^d 2 u — u\u\ 2 + V(x)u , V(x) = W(hx) . 

For V = this equation is completely integrable - see for instance [3j. One of the most 
striking consequences of that is the existence of exact TV-soliton solutions: 

u(x, t) = q N (x, a + tv,v,9 + |(/x 2 + v 2 ), /i) , 

a,v e R N , 9 G (R/2nZ) N , 
where the construction of = 9at(x, a, v, 9, /i) will be recalled in §[2j 
When V ^ and 

u(x, 0) = ^at(x, a, v, 9, /i) , 
the exact dynamics (1.7) is replaced by 

(1.8) u(x, t) = q N {x, a{t),v{t), 0(t),fj,(t)) + 0(h 2 ) , 

where the precise meaning of the error, and its optimality, will be described below. The 
parameters of the multisoliton approximation solve the system of ordinary differential equa- 
tions: 

Vj = -^ l (d a V N + Vjd 0j V N ) , cbj = Vj + ^ x d v V N , 
fij = de 3 V N , ^ = v 2 /2 + Mj 2 /2 + ^\d Vj V N - d N V N , 

and where 

def 1 /* 

Vn(cl,v,6, [jl) = - / ]/(x)|gAr(^, tt, 'U, 0, /x)| 2 rfx . 

^ ./TCP 



Although somewhat complicated looking, the equations ( |1.9[ ) have a natural interpretation 
in terms of Hamiltonian systems: they are the Hamilton- Jacobi equations for the full 
Hamiltonian of (1.6) restricted to the symplectic 4iV-dimensional manifold of TV-solitons - 
see ^3] for details. Of course when V = the solutions correspond to the exact solutions of 
fll.7p . The mathematical results of [7], [8] suggest that the approximation (1.8) is valid up 
to a (rescaled) Ehrenfest time (|1.4[): 



Q holds for < t < C\og{l/h)/h. 

In other words, the equations ( |1.9| ) give the minimal exact effective dynamics valid up to 
the Ehrenfest time log(l//i)//i, where h is the parameter controlling the small variation of 
the potential, see (3.1). In this work, we show that the approximation errors 0(h 2 ) and the 
Ehrenfest time bound are sharp. See [TT] for a survey of soliton dynamics under integrable 
systems that have been perturbed. 

One motivation for this study is the experimental and theoretical investigation of soliton 
trains in Bose-Einstein condensates [14J. We show that the effective dynamics described in 
§[3] is in qualitative agreement with the behaviour of the matter- wave soliton trains - see 
Figure [2| 

The paper is organized as follows: in ^2] we recall the construction of TV-soliton solutions 
for V = and in ^3| the Hamiltonian structure of the equation and the derivation of the 
effective equations of motion. In ^4] we compare the effective dynamics to the behaviour 
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Figure 2. Four solitons with alternating phases bunching up and spreading 
out in the potential V(x) = (x/2) 6 . The full solution is plotted on the right 
with a bird's eye view. The figures on the left are snapshots of that solution. 
Due to their alternating phases, the solitons repel and never pass through 
each other. 



of solutions to (1.6) and draw some quantitative conclusions. Specific potentials similar 
to those in [14J are then discussed in §[5| We investigate effective dynamics for the mKdV 
equation in Finally, in ^7] we describe the numerical methods and compare to other 
possible approaches. 



2. 7V-SOLITONS FOR CUBIC NLS 



When V = 0, we recover the nonlinear cubic one dimensional Schrodinger equation, 
which has TV-soliton solutions with explicit formulas that we now recall - see [3] for a 
detailed presentation of this completely integrable equation. 

We will construct functions ^(x) that depend on 4N parameters: positions, velocities, 
phases, and masses: 

(2.1) q N (x) = q N (x,a,v,9, fi), a,v,9eR N , fi<E(0,oo) N . 
Put 

(2.2) Xj = Vj + ifij , 7j (x) = e^e^-^e^ , 
and define matrices 

M{x) e R NxN , Mi(x) G R^+Dx^+D , 
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-100 exp(cos(x)) 




1000 



-100 sech(5x) 2 + x 4 



300 cos(xr 





Figure 3. A gallery of potentials used for the numerical experiments. Since 
the solitons in the experiments have width approximately 1/10, the inter- 
esting potentials should have size approximately 100. This is suggested by 



the rescaling (3.9). The potentials vary on a scale comparable to 1, hence 



the effective h is approximately 1/10. The exception is the upper right plot 
where we intentionally chose a potential which will exhibit some failures of 
effective dynamics. In the analysis of errors, for instance in Figure |6j only 
relative sizes of h matter. 



by 

(2.3) 

where 

(2.4) 

Finally, 

(2.5) 



M jk (x) 



1 + 7 i (x)7fc(x) 
Xj — Afc 



M l 



M(x) 7 
1 



7 = [7i>" - ,7jv] T , 1 = [!,••• ,1] € 



def detM^x) 
q N (x) = 



det M(x) 



Remarkably, this gives a solution to (1.5) with V = 0, the iV-soliton solution: 

t 



(2.6) 



As one can see from the formula, some restrictions on the parameters apply, see [3]. 
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3. Effective Dynamics Equations 



We consider potentials defined on K. that are slowly varying in the sense that 

(3.1) V{x) = W{hx) 
where W(x) is C 2 in x ) and 

\d k x W(x)\ < C(l + \x\) N ,k < 2. 

where C and TV are independent of h. This means that h is the parameter controlling the 
slow variation of V. 

To obtain an effective dynamics for the evolution we use the Hamiltonian structure of 
the equation. In the physics literature an approach using Lagrangians is more common - 
see for instance Goodman-Holmes- Weinstein [5j and Strecker et al [14]. In the mathematics 
treatments [4j,[8j,[7j the Hamiltonian approach was found easier to use, which we follow 
here. 

The basic claim is that an approximate evolution of is obtained by restricting the 
Hamiltonian flow generated by the Gross-Pitaevskii equation to the manifold of TV-solitons 
described in ^2} The Hamiltonian associated with the Gross-Pitaevskii equation is 

(3.2) H v (u) = j {\d x u\ 2 ~ M 4 ) dx+^Jv\u\ 2 
with respect to the symplectic form 

(3.3) uj{u 1 v) = Im J uv . 

The manifold of solitons, Mjv, is 4j 7 V-dimensional and equipped with the restricted sym- 
plectic form given by the sum of forms for single solitons: 

N 

(3.4) ujm = w\m = ^^(fijdvj A dctj + Vjdfij A daj + d9j A dfij) . 

3=1 

Hy restricted to Mn is 

(3.5) H N = f H v \ MN (a,v,0,n) = ^ f - + V N (a,v,6, /i) , 

j=i \ / 



(3.6) 



where Vn(cl,v,9, ji) 



def 



V(x)\qN(x, a, v, 9, /x)| 2 dx . 



The effective dynamics is given by the flow of the Hamilton vector field of Hn on the 
manifold M N . That vector field, S^, is defined using the symplectic form (3.4): 



(3.7) 



dH N — ujm(', ^h n ) • 
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A computation based on this gives the following ordinary differential equation for the 
parameters a,v,9 and /i, called the effective dynamics: 



(3.i 



Vj = -fjL^d^Htf - /i^Vjde^N = -n]\d aj V N + Vjd 0j V N ) , 
cij = n~ l d Vj H N = vj + \i~ x d V) V N , 
fij — d0jH N = dejVjsi , 

Qj- = ii~ l Vjd Vj H N - d H H N = v]/2 + /j?/2 + ^~ l Vjd Vj V N - d H V N . 



We remark that one can scale the Gross-Pitaevskii equation (1.6) in the following way: 
Consider any function u(x, t), scaling parameter a, let x — ax , t — a 2 t, and define the new 
function 



(3.9) 



u(x, t) = f — u(x, t) . 
a 



Then if u(x,t) satisfies (1.6) with the potential V(x), u(x,t) also satisfies (1.6) with the 
new potential 



(3.10) 



V(x) 



a 2 \a 



This means that if we deal with a soliton of width comparable with a, the potentials for 
which interesting dynamics should appear should have size approximately a~ 2 and the 
slowly varying factor replaced by h/a. 



The effective dynamics equations (3.8) scale similarly: if (a(t),v(t),0(t), fi(t) ) sa tisfies 



(3.8) and we define x, £, and V as above, then ( a(t), 0(t), jl(t) ) also satisfies (3.8) with 



(3.11) 



a(t) = aa(t) , v(t) = — , 9(t) = 0{t) , fl(t) = 



a a 

The scalings (3.9) and ( 3.11[ ) are related in the following way: if 

u(x, t) = q N (x, a(t),v(t),0(t),fj,(t)) , 

then 



4. Comparison of effective and exact dynamics 



For given values a ,^o^o,Mo i n ^ N \ we consider the solution u(-,t) of (1.6) with ini- 
tial data <7jv(*j a , Vo,0 , /x ) and the solutions a(t),v(t),0(t), of the effective dynamics 
equations ( |3.8[ ) with initial values a , ^o, Mo- In the following discussions we will refer to 
u(-,t) as the exact solution and q^(- , a(t) , , 0{t) , /x(t)) as the effective dynamics. 

Holmer and Zworski [8] proved that in the case N — 1, 
(4.1) |K, - a(f), 0(f), = ^ 2 "' . for * < ^glM , 



Exact solution; N = 1 



Exact solution; N = 3 
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Effective Dynamics; N = 1 Exact solution; N = 2 



Effective Dynamics; N = 2 




Effective Dynamics; N = 3 



Exact solution; N = 4 



Effective Dynamics; N = 4 




Figure 4. A side- by-side comparison of the absolute value of the exact 



solution of (1.6) versus the effective dynamics (3.8) for 1,2,3, and 4 solitons 
with potential W(x) = —100 sech 2 (5x) + 10x 4 . The sharpness of the sech 2 (5x) 
term creates clearly visible discrepancy between the two solutions. 



where 5 G (0, 1/2) can be chosen, and where C depends only on the potential and initial 
velocity of the soliton, but not on 5. The H 1 norm measures the size of the function and 
its spatial derivative in L 2 : 



1 2 def 

\H 1 — 



l'117/l = \W\\ 2 L 2 + ||<9^|| 2 2 , 



\ 2 L 2 = f / \v(x)\ 2 dx . 



This norm measures the energy of the solution. 

The limiting time \og(l/h)/h is the Ehrenfest time discussed int ^TJ 

It is expected that this result also holds for TV > 1. This is suggested by [7j, which 
proves the analagous theorem for the modified Korteweg-de Vries (mKdV) equations for 
case N — 1, 2, see ^6] below. However, the methods of [7j do not fully apply to the case of 
the Gross-Pitaevsky equation (1.5). Also, even in the case of mKdV and = 2, multiple 
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soliton interactions are not theoretically understood. All these considerations provided a 
strong motivation for this numerical study. 



We note that for N > 1, it was conjectured in [7J that the error bounds (4.1) will hold 
not only in the H 1 norm, but for the H N norm, which measures the the size of a function 
and its first TV derivatives, where TV is the number of solitons. This has been proven in the 
mKdV case with TV = 2 [7j. For our numerical experiments, we consider only the H 1 norm. 



We present numerical simulations to show that the result (4.1) holds for TV > 1 in the 



following three sections: In §4.1| we choose initial data and two potentials that demonstrate 
the power and limitations of the effective dynamics equations, regardless of the number 

we 



of solitons. Using the same initial data and one of the potentials from §4.1[ in §4.2 
verify that the 0(h 2 ~ s ) error estimate in (4.1) holds as h — >> for a fixed time interval. We 



then turn to the \og{l/h)/h timescale, or Ehrenfest timescale, in £4.3 to show that it is the 



appropriate timescale for which we can expect (4.1) to hold. 



4.1. A numerical case study. We consider initial data gjv(-, a/v, v^, 9n, Mat), where = 
(ai, . . . , aw) , A = 1, 2, 3, 4 and v^, 9n-> and JIn are similarly defined with 

(ai, a 2j a 3 , a 4 ) = (-1, -1.5, 0, 1) 
Oi, v 2) v s , v 4 ) = (-2, 0, 3, 0) 
( " } (^,^2,^3^4) = (vr/3, 0,-3, -5) 

(/ii, /i 2 , /x 3 , fjb 4 ) = (17, 25, 23, 19) 

The positions and masses are chosen to satisfy a numerical requirement that ( • , cljsi , vn -> @n 
is close to outside of (— 7r, 7t), our numerical domain (see Rescaling the solution as in 



(3.9) or enlarging the numerical domain allows for data that does not satisfy this numerical 



requirement. The initial data is otherwise chosen arbitrarily. 
We first consider the potential 

Vi(x) = -100e cosx , 

see Figure [3| The factor —100 is chosen to create a deep enough well so that the solutions 
remain in (— 7r,7r), but the potential is otherwise chosen arbitrarily. 

We compute the exact solution and the effective dynamics solution for = 1,2,3,4 
up to time t — 1, which is chosen to allow for multiple soliton interactions. We plot the 
solution for 4 solitons in Figure [TJ where we observe very little difference between the exact 
and effective dynamics solutions. An equally small amount of discrepancy between the 
solutions was observed for N = 1, 2, 3. 

Next, we consider the same initial data as above with the potential 

V 2 (x) = -100sech 2 (5x) + 10 x 4 , 

see Figure [3] This potential is chosen to be outside of the slowly varying regime for which 
the effective dynamics give good approximations. This is due to the — 100sech 2 (5x) term, 
which creates a sharp dip roughly the width of the solitons we are studying; thus we do 
not expect the exact solution to maintain its soliton structure very well. This causes the 



10 



TREVOR POTTER 








; tr| 


rr 


- Exact Solution 

Effective Dynamics 
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X 

t = .658 


12 3' 







-rl 




1 Exact Solution 

Effective Dynamics 



Figure 5. The plot on the left is a different view of the exact solution with 
3 solitons shown in Figure [4| The plot on the right compares the real parts 
of the N = 3 exact solution with the effective dynamics solution at times 
t = 0.35 and t = 0.658. 



clearly visible discrepancies between the effective dynamics and the exact solution in Figure 
[4j The 10x 4 term ensures the solutions remain on the interval (— 7r,7r). 

We compute the solutions for N = 1, 2, 3, 4 up to time t = 0.7, which is again chosen to 
allow for multiple soliton interactions. Figure [4] displays these solutions and demonstrates 
that the effective dynamics captures the true motion, regardless of the number of solitons 
and regardless of multiple soliton interactions. In the experiments presented in this paper, 
we only consider TV < 4, but we have observed good agreement between the exact solution 
and the effective dynamics for N < 7. We did not investigate futher due to increasing 



computational time needed to solve (3.8). 



We note that for TV > 2 the phases of the solitons are crucial in determining the inter- 
action between solitons. In Figure [4j we see that for TV = 3, at approximately t ~ 0.65, 
two solitons that appear to bounce off each other in the exact solution instead appear to 
cross in the effective dynamics. This discrepancy seems to be due to differences between 
exact phases and effective phases. In Figure [5] we are able to see large deviation in the 
phases between the exact solution and effective dynamics by comparing the real part of the 
solutions. 



4.2. Quantitative study of the error as h —> 0. We investigate the 0(h 2 ~ 5 ) error 
between the exact solution and the effective dynamics on a fixed time interval. The estimate 
that gives rise to (4.1) is 

(4.3) t) - q N (; a(t),v(t), 9(t)^(t))\\ H i < Ch 2 e Cht . 

If t = 5log(l/h)/(Ch), then the RHS reduces to Ch 2 ~ s . When dealing with fixed length of 
time or even time of size 0(1/ h) the RHS is 0(h 2 ) and that form of error will be shown to 
be optimal. 
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Convergence of solution of effective dynamics equations to exact solution as h tends to zero 




Figure 6. A log-log plot of the H 1 error, relative to the H 1 norm of the 



initial data, between the exact solution of |1.6| and the TV-soliton evolving 
according to the effective dynamics equations |3.8| as function of h. Here, 
the potential is V(x) = W{hx), where W(x) = — 100sech 2 (5x) + 10x 4 . For 
smaller values of the slope of the lines approaches 2, in agreement with 
the theoretical upper bound on the error in (11.81). 



We reconsider our second potential from §4.1[ but add the slowly varying parameter, h: 

V(x) = W(hx) , W(x) = -100sech 2 (5x) + 10x 4 , 

and explore the H 1 error, relative to the H 1 norm of the initial data, between the exact 
solution and effective dynamics as h —> 0. We expect that as h becomes small enough, the 
equation will enter the slowly varying regime and display 0{h 2 ) error. Indeed, the log-log 
plot in Figure [6] demonstrates the error is bounded by Cjsfh 2 as h —> 0, where the constant 
Cn varies only slightly between different values of N. 

We fit the data from Figure [6] to a line using the 6 smallest values of h. 



N 


1 


2 


3 


4 


Slope 


1.86 


1.76 


1.91 


1.86 


Cn 


-1.07 


-1.46 


-1.27 


-1.37 



Thus, we conclude that the error is approximately 0(h 2 ). 

4.3. Ehrenfest time. We now inv estig ate the length of time for which the effective dynam- 
ics approximation is accurate. In (4.3), we recalled that the error is bounded by Ch 2 e Cht 



and hence the approximation breaks down at the Ehrenfest time, 

T(h)~l 0g (l/h)/h. 

We have already verified for small h and fixed time this error behaves as 0(h 2 ). Thus we 
focus on observing exponential growth in the error as a function of time, and verifying that 
it is of the form 0(e Cht ). 

For this we must choose a potential and initial data to exhibit exponential instability. 
We are motivated by Newton's equations for V(a) = —a 2 /2: 

(4.4) a — v , v — a , v = v cosh t + a sinh t , a = a cosh t + v sinh t . 
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error between exact solution and effective dynamics 



linear dependence of B on h 




co 40 



10 



0.15 0.2 
time 



0.6 0.6 
h 



Figure 7. The plot on the left shows H 1 error, relative to the H 1 norm 
of the initial data, between the exact solution of (1.6) and the effective 



dynamics for a single soliton sliding down the concave potential V(x) = 
W{hx) , where W(x) = 1000 x 2 . The error is plotted as a function of time 
and for several values of the parameter h. On the right sight, B is plotted as 
a function of the parameter h ) when the errors from the plot on the left are 
fitted to a curve of the form A(e Bt + C). We expect B to depend linearly on 
h. 



In this case, we have exponential instability of classical dynamics. This suggests choosing 
potentials with a non-degenerate maximum and working near the unstable equilibrium 
points. 



Hence we will investigate solutions to (1.6) with potential 



V(x) = W(hx) , where W(x) = -1000 x 2 



Figure [7] below demonstrates exponential divergence between the exact solution to (1.6) 



and the effective dynamics for several values of h and a single soliton initial condition 
qi (x,. 1,0, 0,15). 

We fit the plots shown in Figure [7] to a function of the form A(e Bt + C), for the time 
period when the soliton's position was between x = .15 /h and x = 1.2/ h. This range was 
observed to be a region where exponential increase dominated the error and before the 
soliton approached the numerical boundary. 



In Figure [7] we observe a linear dependence of B on h : in agreement with (4.3). This 



indicates that for certain potentials the Ehrenfest time C\og(l/h)/h is the appropriate 
bound for the length time we expect the effective dynamics to give a good approximation 



to (1.6). We note that in our experiments with other potentials we often observe a linear 
increase in error which would correspond to a timescale of C/h 2 instead of the Ehrenfest 
time C \og(l/ h)/h. 
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Figure 8. Surface plots of the error between the exact solution of ( |1.6D and 
the effective dynamics for a single soliton sliding down the concave portion 
of the potential V(x) = W(hx), where W(x) = 1000 x 2 , as in Figure [7[ 
We have plotted the absolute value of the difference between the spatial 
derivatives between the two solutions. 



5. Application to Bose-Einstein Condensates 

Strecker, Partridge, Truscott and Hulet [14] discovered that Bose-Einstein condensates 
form stable soliton trains while confined to one-dimensional motion. When set into motion 
in a suitably chosen optical trap, a Bose-Einstein condensate forms multiple soliton forma- 
tions which exist for multiple oscillatory cycles without being destroyed by dispersion or 
diffraction. We can observe this same behavior numerically, using the effective dynamics 
equations. We choose a potential of the type described in [2] 

= (i) 6 

and set TV = 4, which was the most frequent case in their experiment. Strecker et al inferred 
that the repulsive behavior of the solitons indicated alternating phases. Their argument 
was based on considering a certain reduced Lagrangian. 

In our numerical experiment we put 9 = (0, 7r, 0, tt) and then set the four solitons in 
motion with the same velocity near the center of the potential. Similar to [2], we observe 
bunching and spreading of the soliton train for several oscillations. See Figure [2} 

6. Effective dynamics for the mKdV equation 



The mKdV equation 
(6.1) d t u = -d x (d 2 x u + 2u 3 ), 
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-2 




Figure 9. The left plot shows a side- by-side comparison of the exact solu- 
tion of the mKdV equation ( |6.2[ ) and the effective dynamics solution for 3 
solitons with potential b(x) = 300 cos 2 x. No discrepancy between the two 
solutions is visible. The right plot displays the exact solution from a different 
angle. 



like the nonlinear Schrodinger equation, has soliton solutions and a Hamiltonian structure. 
Holmer, Perelman, and Zworski [7J derived effective dynamics equations for the mKdV 
equation with a slowly varying potential 



(6.2) 



d t u = —d x {p 2 x u — b(x 1 t)u + 2u 3 ), t) = b (hx, ht) 



and proved a result analogous to the (|4.1|) for TV = 1,2: the H N error between the solution 



of (|6.2|) and its associated effective dynamics with TV-soliton initial data is bounded by 

1 



(6.3) 



Ch z e 



2 CM 



c 

for t < — log , 
h h 



Similarly to ^42, we have conducted a numerical study verifying that the H 1 error is 
0(h 2 ) as h for multiple soliton initial conditions. See Figures ^ and 10, 



7. Numerical methods 



We now describe the numerical methods we employ to compute the Gross-Pitaevskii PDE 
(1.6) and the ODE (3.8) arising from the effective dynamics. When comparing a solution 



of (1.6) with a solution of (3.8), we refine our numerical solutions until the error between 



sucessive refinements of solutions to the same equation is several orders of magnitude 
smaller than the error between solutions of the two equations. 



EFFECTIVE DYNAMICS FOR iV-SOLITONS 



15 



Convergence of solution of effective dynamics equations to exact solution as h tends to zero 




Iog10(h) 



Figure 10. A log-log plot of the H 1 error, relative to the H 1 norm of the 



initial data, between the exact solution to the mKdV equation (6.2) and the 



AT-soliton evolving according to the effective dynamics, as a function of h. 
For smaller values of h, the slope of the lines approaches 2, in agreement with 
the theoretical upper bound on the error. The theoretical upper bound has 
only been proven for TV = 1,2, but this figure gives evidence that it holds for 
all N. 



Numerically solving the ODE arising from the effective dynamics (3.8) necessitates com 



puting qjsr(x,a,v,9, //) and its derivatives with respect to the parameters a, v, #,and \i effi- 
ciently. For this we note that an equivalent definition of q^ in (2.5) 



(7.1) 



-1M~ 



7: 



where 1, M, and 7 are as in (2.5). Since iM is Hermitian, M lr y can be efficiently computed 
using the Cholesky factorization. Differentiating q^ numerically, for larger values of TV, is 



too costly. Instead we used (7.1) to obtain explicit formulas for the derivatives of q^ and 



again used Cholesky factorizations to efficiently compute them. With this we compute the 
integrand in (3.6), and then numerically integrate it using the trapezoidal method. Once 



we can efficiently compute the RHS of the effective dynamics equations (3.8), the standard 



fourth-order Runge Kutta method was found to be suitable to solve to the ODE. 



In order to solve the PDE (1.6) we used a Fourier spectral method to study the evolution 
on the numerical domain (— 7r, tt). This requires our solution u{x, t) to be periodic in space, 
so we choose initial conditions such that u{x, t) decays to zero, to machine precision, before 
the endpoints — tt and tt. Arbitrary initial data can be handled by either extending the 
numerical domain or rescaling the equation (see ( |3.9D ). 

One difficulty arises in that a non-trivial potential W(hx) cannot be periodic for all h. 
However, the potential W{hx) need not be periodic on (— 7r,7r) so long as the product 
W{hx)u{x ) t) is periodic on (— 7r,7r), which is achieved if u(x,t) decays fast enough at 
the endpoints — tt and n. If W(x) is periodic on (— 7r,7r) and one wishes to consider a 
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solution u(x,t) that doesn't decay before the endpoints — tt and 7r, the rescaling (3.9) may 
be employed with a = h: 



(7.2) 



ti(x, t) = f \u{x, t) 

h 



Then if u(x,t) satisfies (1.6) with periodic potential V(x) 
(fL6l) with potential V(x) = W(hx). 



W(x)/h 2 , u(x,t) also satisfies 



This rescaling also makes it clear that as h — >> 0, a soliton solution becomes sharper 
relative to the potential. This requires higher resolution in order to apply our numerical 



method to solve the PDE (1.6), while the effective dynamics equations (3.8) are unaf- 



fected. Indeed, our numerical experiments confirmed that increased computational effort 
was needed to resolve the PDE as h 0, but not the effective dynamics ODE. 



We now describe the method to solve a general solution of (1.6) on a periodic 

domain with periodic initial data and potential V(x). The Fourier modes iik(t) of a solution 
u(x,t) to (1.6) evolve according to 



(7.3) 



d t Uk = --k 2 u k + i{uv) k , 



V \u\ 



V 



Discretizing space and replacing the Fourier Transform with the Discrete Fourier Transform 
gives rise to a finite dimensional system of ODE, which we now represent in the general 
form 



(7.4) 



u t 



Cu + M{ 



u 



where C is a stiff linear transformation corresponding to the first term of ( |7.3[ ) (represented 
by a diagonal matrix in our case) and N is a non-linear operator from the second term 
of ( |7.3[ ). To solve (7.4) we compared the fourth order implicit-explicit (IMEX) method 
ARK4(3)6L[2]SA proposed by Kennedy and Carpenter [10] with the exponential time dif- 
ferencing (ETD) method ETDRK4 proposed by Kassam and Trefethen [9]. The IMEX 
scheme update formula is 

(7.5) u n+1 = u n + A*(&i(fci + h) + • • • + b s (k s + l s )) , 
where At is the time step and ki and U are chosen such that 

(7.6) k % = C{u n + A*(^ifci + • • • + A is k s + Anh + ••• + A is l s )) 



(7.7) 



h = M{u n + AtiAah + • • • + Aisk 8 + Anh + • • • + A is l s )) 



Here A, A are s x s lower triangular matrices with A having zeros along its diagonal.. This 
allows us to solve for the ki and li one stage at a time, only inverting the diagonal linear 
operators / — AtAuC. The implicit treatment of the C term mitigates the stiffness arising 



from the k 2 factor in (7.3), while the li can be computed explicitly, so non-linear equations 



involving J\f need not be solved. 

The ETD method, on the other hand, uses an exact formula for obtaining the next step 
u n +i from u n based on solving the linear portion exactly: 

"At 



(7.8) 



u n+1 e 



CAt 



/ 

Jo 



e- LT N{u{t n + T),t n + r) dr 
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Convergence of numerical solution of PDE as the timestep tends to zero 



Convergence of numerical solution of PDE as function of computational time 





ETDRK4 
ARK4(3)6L[2]SA 



log 10 (timestep) 



Iog10(computational time in seconds) 



Figure 11. Log-log plots of the convergence of the fourth order schemes 
ETDRK4 and ARK4(3)6L[2]SA as a function of the timesteps and compu- 
tational time, respectively. ARK4(3)6L[2]SA is slightly more efficient per 
timestep, but ETDRK4 is significantly more computationally efficient. 



The integral in (7.8) can then be numerically approximated using matrix exponents of C 
and evaluations of N. Thus as with the IMEX method, we do not need to solve non-linear 
equations, and computations involving C (namely computing e AtC ) are efficient because 
C is diagonal. Stiffness is mitigated by solving the linear portion of (7.4) exactly. We 



found that the ETDRK4 scheme computed a solution of a desired accuracy nearly twice 
as fast as the ARK4(3)6L[2]SA scheme. While the ARK4(3)6L[2]SA scheme had a slightly 
smaller error rate per step, more computations per step made it significantly less efficient. 
Neither method demonstrated any instability in the range of step sizes required for our 
solutions. Below, we plot the convergence of the two schemes as the timestep goes to zero. 



To obtain the results in the figures, we used the same potential and initial data as in §4.1 
W(x) = -100sech 2 (5x) + 10x 4 . 
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